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INTRODUCTION 


Current  world  conditions  require  that  the  Navy  has  the  ability  to  rapidly  make 
geoacoustic  surveys  in  areas  that  were  previously  of  little  interest.  Past  and  current 
survey  methods  are  cumbersome  and  slow.  This  report  addressed  this  issue  by 
developing  a  rapid  method  of  solving  for  the  environmental  properties  within  a  survey 
area.  This  method  can  be  applied  to  a  variety  of  geoacoustic  data  such  as  multichannel 
marine  seismic  data  but,  in  particular,  we  wish  to  acquire  geoacoustic  data  with  an  array 
of  sonobuoys  and  SUS  charges  that  are  fairly  randomly  arranged  (i.e.,  no  rigorous 
geometric  constraints  required).  Then,  using  these  data,  we  solve  for  the  subbottom 
geoacoustic  parameters  such  as  compressional  wave  velocity  (Vp)  illustrated  in  Fig.  1. 
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Fig.  1 — Cross  section  of  subbottom  profiling  experiment  with  floating 
receivers  and  multiple  pulse  sources.  Vp  values  are  plotted  at  the  right. 

To  find  a  solution  by  our  method  we  need  to  search  over  a  large  model  space  (Fig.  2). 
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Fig.  2 — Illustration  of  the  search  over  model  space  to  find  the  true  values 
of  Vp.  The  left  figure  illustrates  the  true  but  unknown  values  of  Vp^  and 
the  right  figure  illustrates  the  lower  and  upper  bounds  of  the  Vp  and  layer 
thickness  values  for  the  search. 

Using  a  relatively  simple  6-layer  model  for  continental  shelf  sediments  with 
resolutions  possible  with  high-frequency  multichannel  geoacoustic  systems  (1  to  10  m  for 
layer  thickness  and  10  to  100  m/s  for  Vp)  there  are  1.1  *  10^^  solutions  in  this  model 
space.  Sampling  one  solution  every  3.25  s  (the  actual  rate  on  a  Sun  SPARC  10  desktop 
workstation),  it  would  take  34,000,000,000  years  to  do  a  complete  search.  Using 
Simulated  Annealing  (SA),  a  directed  search  method,  we  found  a  solution  with  maximum 
errors  of  less  than  5%  and  average  errors  of  1%  to  2%  for  both  layer  thicknesses  and 
sound  velocities  in  2  h  and  25  min  (Fig.  3).  If  this  inversion  was  run  on  a  massively 
parallel  computer,  then  we  would  expect  run  times  of  a  few  minutes. 

SA  is  a  nonlinear  inversion  method  (see  Appendix)  that  searches  model  space  in  a 
pseudo  random  fashion  that  is  directed  by  a  statistical  selection  criterion  (Kirkpatrick  et 
al.  1983).  The  search  over  model  space  is  done  by  altering  one  variable  in  the  current 
model  at  a  time,  creating  a  new  model,  and  then  calculating  synthetic  data  with  a  forward 
model  routine.  Synthetic  data  are  compared  to  the  field  data  using  an  objective  function. 
A  variety  of  objective  functions  can  be  used,  such  as  cross  correlation  or  residual  (Sen 
and  Stoffa  1991).  If  the  new  model  produces  a  lower  objective  function,  then  this  new 
model  is  chosen  as  the  current  model  and  the  process  is  repeated.  If  the  new  model 
produces  a  higher  objective  function,  then  it  is  occasionally  chosen  as  the  current  model 
by  a  method  that  resembles  the  statistical  energy  distribution  of  atoms  in  a  solid.  At  a 
high  "temperature,"  the  chance  of  choosing  the  solution  with  a  higher  cost  function  is 
much  higher  than  at  a  low  "temperature."  This  process  is  begun  at  a  high  temperature  and 
repeated  hundreds  or  thousands  of  times  as  the  temperature  is  slowly  reduced  in  a  process 
that  is  mathematically  analogous  to  annealing  a  solid  object.  If  the  temperature  is 
reduced  slowly  enough  over  the  proper  range,  then  the  solution,  or  object,  will  "anneal" 
into  a  single  crystal,  which  is  the  lowest  energy  state  possible,  and  the  synthetic  data 
calculated  will  be  the  best  possible  fit  to  the  field  data  (van  Laarhoven  and  Aarts  1987). 
See  the  appendix  for  more  details  about  SA. 
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Current  Status  of  the  Inversion  Method 


The  working  SA  algorithm  is  based  on  the  classical  Metropolis  algorithm  (Metropolis 
et  al.  1953)  (see  Appendix)  but  has  significant  differences  in  all  major  steps.  Two  values 
used  throughout  the  algorithm  are  the  objective  function  and  the  temperature.  The 
objective  function  is  used  to  measure  the  goodness-of-fit  between  the  synthetic  data 
calculated  from  the  test  model  m  and  the  field  data  (with  a  small  value  being  a  better  fit 
than  a  large  value),  and  the  temperature  is  a  parameter  used  in  selecting  or  rejecting  a  bad 
fit  and  to  adjust  the  size  of  the  model  parameter  search  window.  We  can  use  any  of  five 
different  objective  functions;  the  zero  lag  cross  correlation,  a  modified  cross  correlation 
that  retains  amplitude  information,  the  residual  (sum  of  the  differences),  the  residual 
squared,  and  the  square  root  of  the  residual. 

We  use  a  temperature  dependent  search  distribution  and  a  layer  dependent 
temperature.  A  temperature  dependent  search  distribution  gives  a  much  finer  search 
pattern  near  the  minimum.  The  layer  dependent  temperature  fits  the  top  layers  first  and  is 
more  efficient  for  a  model  with  many  layers  as  it  results  in  a  linear  increase  in  the  run 
time  relative  to  the  number  of  layers  rather  than  a  geometric  or  exponential  increase  in 
time. 

We  tested  the  SA  algorithm  by  inverting  synthetic  data.  The  synthetic  data  are 
simplified  and  noise-free  examples  of  what  is  expected  from  real  geological  sites.  The 
synthetic  data  simulate  field  data  collected  by  sonobuoys  using  SUS  charge  sources  as 
well  as  field  data  from  high-resolution  multichannel  systems  such  as  the  Deep  Towed 
Acoustics-Geophysics  System  (Gettrust  and  Ross  1990).  We  can  invert  seismic  or 
geoacoustic  data  in  the  tau-p  domain,  the  t-x  domain  (time  and  offset,  the  most  direct 
observation),  or  the  to-x  domain  (frequency  and  offset)  with  any  of  five  different 
objective  functions.  Geoacoustic  data  are  normally  collected  in  the  t-x  domain  and 
transformed  to  the  tau-p  or  the  co-x  domains  if  desired.  However,  a  reliable  tau-p 
transformation  requires  more  offsets  (i.e.,  sonobuoy-SUS  pairs)  than  necessary  for  the 
inversion  done  in  the  t-x  domain. 

Figures  3  and  4  are  examples  of  simulated  field  data  and  inversion  solutions.  The 
example  shown  in  Fig.  3  uses  the  expected  geoacoustic  profile  for  a  thick  sedimentary 
basin  such  as  a  continental  shelf,  specifically  the  Nile  Fan  in  the  Mediterranean  Sea. 
These  synthetic  data  are  shown  in  the  tau-p  domain  (intercept  time  and  ray  parameter,  a 
plane-wave  decomposition  of  the  wave  field),  which  makes  for  a  faster  inversion 
calculation.  Figure  3a  shows  the  synthetic  "data"  in  tau-p  space.  Figure  3b  shows  the 
residual  or  the  difference  between  the  "data"  and  synthetics  calculated  from  the  solution 
or  fit  found  by  the  inversion  calculation.  The  low  amplitudes  here  indicate  that  the 
solution  model  is  close  to  the  true  model  (the  fit  is  good).  Figure  3c  shows  the  cooling 
curve  or  the  progress  of  the  fit  as  the  calculation  proceeds  from  a  high  temperature  and 
large  objective  function  on  the  right  side  to  a  low  temperature  and  small  objective 
function  on  the  left  side.  The  objective  function  is  best  defined  so  that  a  perfect  fit  gives 
a  value  of  0  and  a  completely  random  guess  gives  an  average  value  of  about  1.  Figure  3d 
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compares  the  model  used  to  generate  the  "data"  with  the  solution  found  by  inversion.  The 
solid  line  profile  was  used  to  calculate  the  "data"  and  the  dashed  line  profile  is  the 
solution  from  the  19th  run  of  our  SA  inversion  method.  We  solved  for  layer  thickness 
and  Vp.  The  closeness  of  the  two  profiles  demonstrates  that  this  solution  is  good. 

Figure  4  shows  an  interesting  geological  situation  found  with  a  multichannel  seismic 
marine  survey;  that  of  one  or  two  thin  layers  having  high  over  otherwise  normal 
sediment.  Figure  4a  shows  the  synthetic  "data"  in  t-x  space.  The  offsets  here  are 
regularly  spaced  but  need  not  be.  Irregular  spacing  as  from  sonobuoys  and  SUS  charges 
will  work  well  in  this  domain.  Figure  4c  shows  the  cooling  curve  that  shows  a  decreasing 
objective  function  at  the  lowest  temperature  suggesting  that  an  even  better  fit  could  easily 
be  achieved.  However,  the  fit  is  very  close  and  the  residual  is  very  small  (Fig.  4b). 
Figure  4d  compares  the  model  used  to  generate  the  "data"  with  the  solution  found  by 
inversion.  The  resolution  and  frequency  are  much  higher  than  in  Fig.  3  so  that  the  thin 
layers  can  be  resolved.  The  only  variable  searched  over  in  the  inversion  was  the  shear 
velocity  (Vg)  giving  only  three  variables.  This  inversion  was  done  without  layer  stripping 
but  with  a  temperature  dependent  search  window.  The  inversion  done  in  the  t-x  domain 
is  very  slow  even  though  the  search  was  only  over  three  variables. 

The  advantage  that  SA  (or  other  nonlinear  inversion  methods)  has  over  linear 
methods  is  demonstrated  in  Fig.  5,  which  is  a  gray  tone  contour  plot  of  the  objective 
function  for  Vj  of  the  first  and  second  layers  for  the  model  shown  in  Fig.  4.  The 
objective  function  used  here  is  the  negative  of  the  cross  correlation  so  that  a  perfect  fit 
has  a  value  of  0  while  random  models  give  an  average  value  of  1.  The  black  area  is 
where  the  objective  function  has  a  value  near  0  and  the  true  Vg  values  are  500  m/s  for 
layer  1  and  650  m/s  for  layer  2.  The  gray  areas  below  it  are  local  minima.  The  highest 
contour  level  on  this  plot  is  at  0.1,  which  is  a  fairly  good  fit  but  this  covers  much  of  the 
model  space.  The  complex  shape  and  multiple  minima  are  a  challenge  for  the  inversion 
method.  A  linear  inversion  method  will  easily  fall  into  one  of  the  local  minima  and  the 
user  may  be  satisfied  with  an  objective  function  of  0.06  to  0.08,  but  will  end  up  with  the 
wrong  solution.  Our  SA  method  easily  found  the  global  minimum  every  time,  even  when 
the  search  was  done  very  rapidly.  Other  objective  functions  calculated  for  the  same 
model  give  similar  plots  with  multiple  local  minima.  Any  model  with  thin  layers  will  be 
nonlinear  and  have  local  minima  in  their  objective  function  meaning  that  a  nonlinear 
search  method  such  as  our  SA  algorithm  is  needed  to  automatically  find  the  optimum 
environmental  model. 


Plans  for  the  Future 

Our  SA  methods  work  well  with  synthetic  data  in  the  tau-p,  t-x,  or  co-x  domain. 
While  the  tau-p  domain  inversion  is  much  faster  than  the  t-x  domain  inversion,  real  data 
are  normally  collected  in  the  time-offset  (t-x)  domain  and  then  transformed  to  tau-p  or 
another  data  domain  for  the  most  efficient  analysis.  A  reliable  transformation  depends  on 
adequate  coverage  of  offsets,  which  is  an  important  experimental  design  parameter. 
Realistic  experiment  designs  never  have  adequate  offset  coverage  for  a  good  tau-p 
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transformation  making  the  subsequent  inversion  in  that  domain  unreliable.  Inversion  in 
the  frequency-offset  (to-x)  domain  avoids  this  offset  coverage  problem  but  is  currently  no 
faster  than  the  t-x  domain  inversion.  We  will  develop  a  much  more  efficient  (O-x  forward 
model  during  the  current  year. 

SA  inversion  produces  a  single  "best"  solution  which,  if  done  properly,  will  be  very 
close  to  the  "true"  solution.  It  will  not  tell  how  sensitive  the  data  are  or  the  solution  is  to 
the  different  parts  of  the  physical  model.  For  example,  typical  sea-floor  data  will  be  very 
sensitive  to  compressional  velocity  differences  in  thick  layers  but  not  in  thin  layers.  This 
information  is  not  included  in  the  solution,  and  the  "best"  solution  may  be  substantially 
different  from  the  "true"  solution  for  the  model  parts  for  which  the  data  are  insensitive. 
The  "Heat  Bath"  method  is  similar  to  SA  but  produces  a  probability  density  function  that 
gives  a  "best  solution"  at  its  peak  but  also  shows  how  reliable  the  solution  is  by  its  width 
(Basu  and  Frazer  1990).  During  the  current  year,  we  will  incorporate  the  Heat  Bath 
method  as  an  option  to  the  inversion  in  all  of  the  data  domains  mentioned  above. 

Many  improvements  will  be  made  to  the  algorithm,  such  as  cleaning  up  the  inversion 
code,  improving  the  data  windows,  making  the  model  search  window  a  function  of  the 
value  of  the  objective  function  and  the  temperature,  using  a  variable  dependent 
temperature,  using  a  better  normalization  of  the  objective  functions,  making  layer 
acoustic  travel  times  independent  of  the  layer  thickness,  using  improved  cooling 
schedules,  and  investigating  the  sensitivity  to  noise.  We  will  investigate  these  and  other 
ideas  this  year  and  incorporate  all  of  those  that  prove  worthwhile  into  our  method. 

This  year  we  plan  to  apply  this  method  on  field  data  such  as  those  collected  during  the 
Critical  Sea  Test-8  experiment  last  spring  and  during  upcoming  cruises  gathering  data 
that  will  be  appropriate  for  our  analysis  methods. 

Summary 


We  developed  and  tested  an  SA  algorithm  for  inverting  synthetic  geoacoustic  data 
from  a  marine  environment.  We  use  several  significant  improvements  in  the  classical 
Metropolis  SA  algorithm  that  improve  the  speed  and  accuracy  of  the  method,  such  as 
using  variable  dependent  temperatures  and  a  temperature  dependent  model  search 
window.  This  method  works  with  three  data  domains  that  are  either  optimal  for  data 
collection  or  forward  model  calculation,  but  not  both.  We  demonstrated  the  effectiveness 
of  the  SA  method  with  a  realistic  but  unexpected  geological  case  that  is  a  nonlinear 
problem  with  several  traps  that  foil  other  inversion  methods. 
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temperature  vp  (m/s) 

Fig.  3 — a.  Synthetic  tau-p  data  from  the  expected  environment  in  the  Nile 
fan  region,  b.  The  residual  of  the  fit  for  the  inversion  solution  from  the 
19th  run.  c.  The  cooling  curve  shows  the  progress  of  the  calculation  from 
high  temperatures  and  large  objective  functions  on  the  right  to  low 
temperatures  and  small  objective  functions  on  the  left,  d  compares  the 
model  used  to  generate  the  "data"  with  the  solution  found  by  inversion. 
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temperature  Yg  (^n/s) 


Fig.  4 — a.  Synthetic  data  in  the  t-x  domain  from  a  region  with  a  high  Vg 
at  the  sea  floor,  b.  The  residual  of  the  fit  for  the  inversion  solution  which 
solved  only  for  Vg.  c.  The  cooling  curve  shows  the  progress  of  the 
calculation  from  high  temperatures  and  large  objective  functions  on  the 
right  to  low  temperatures  and  small  objective  functions  on  the  left.  d.  A 
comparison  of  the  known  geophysical  profile  and  the  solution  from  our 
SA  inversion  method. 
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Fig.  5 — A  contour  plot  of  the  objective  function  CCl  for  Vg  of  the  top  two 
sediment  layers  for  the  model  described  in  Fig.  2.  While  the  objective 
function  varies  from  1  to  0,  only  the  range  between  0. 1  and  1  is  shown  in 
order  that  the  details  are  brought  out.  The  complex  shape  and  multiple 
minima  are  a  challenge  for  the  inversion  method.  A  linear  method  failed 
in  several  tries  by  either  falling  into  one  of  the  local  minima  or  not  finding 
a  large  enough  gradient  to  refine  the  initial  solution. 
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Appendix:  Nonlinear  Inversion  and  the  Metropolis  Algorithm 

Inversion/optimization  methods  want  to  find  the  lowest  value  of 

E  =  ^d  —  G(m)||  where  E  is  the  error,  d  is  the  observed  data,  G  is  the  theory,  and 
m  is  the  earth  model. 


m 


•  At  each  step,  conventional  inversionyoptimization  goes  to  a  lower  value  of  E  (m),  so 
the  answer  may  not  be  the  global  minimum. 

•  In  any  step,  simulated  annealing  (SA)  may  go  to  a  higher  value  of  E.  Lower  values  of 
E  are  preferred,  but  not  so  strongly  preferred  as  to  send  SA  into  a  local  minimum 
from  which  it  will  never  escape. 

•  SA  works  by  sampling  from  the  Gibbs  distribution: 


m 


while  slowly  lowering  the  temperature  T  without  actually  knowing  Pj. 
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The  Metropolis  Algorithm 


1. 

2. 

3. 

4. 


Choose  a  random  starting  model  m  and  a  high  T.  Define  E(m)  so  that  l^niax  '  ^min' 
=  1 


Randomly  perturb  m  i ,  the  first  component  of  m,  to  get  m'.  If  E'^  then  accept  ni  as 

(£-£')/ 

the  new  model.  If  E'>E  accept  m’  with  probability  e 


/T 


Visit  the  remaining  components  of  m  updating  in  a  similar  manner.  One  visit  to  each 
component  constitutes  a  "sweep." 


Reduce  T  slightly  and  go  to  2.  Stop  if  m  has  not  changed  in  the  last  100  sweeps. 
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